# Replication materials for 
# Peisker, Hoffmann, Muttarak (2025)
# Climate news mediates extreme weather effects on climate change concern

packages <- c("tidyverse", "data.table", "fixest")
lapply(packages, library, character.only = TRUE)
lapply(packages, citation)

source("R Files/functions.R")

#### load data and set variable names ####
load("Data/eb_climate_weekly.RData")
setnames(eb_climate_weekly, c("nuts"), c("unit"))

# variables to standardize with within SD
in_vars <- c(
  "cc_eu_concern_w", 
  grep("tmean_", names(eb_climate_weekly), value = TRUE),
  grep("emm_", names(eb_climate_weekly), value = TRUE),
  grep("google_", names(eb_climate_weekly), value = TRUE)
)
# do not standardize
id_vars <- c("year", "week", "unit")
other_vars <- c(
  "season", "spring", "summer", "autumn", "winter",
  grep("event", names(eb_climate_weekly), value = TRUE),
  grep("env_dt", names(eb_climate_weekly), value = TRUE),
  grep("_year_factor", names(eb_climate_weekly), value = TRUE),
  "nutsyear", "nuts0year", "regionyear", "region",
  "t", "lat", "lon", "norm_dist"
)
all_other_vars <- c(id_vars, other_vars)

#### get fixed effects residuals and calc sd ####
within_sd <-
  sapply(in_vars, function(i){
    felm_tmp <- feols(
      f_feols(yvar = "norm_dist", xvars = i, fevars = m_fv),
      eb_climate_weekly, demeaned = TRUE)
    sd_tmp <- c(sd(felm_tmp$X_demeaned, na.rm = TRUE))
    return(sd_tmp)
  }) %>%
  data.frame() %>%
  rownames_to_column
names(within_sd) <- c("variable", "within_sd")
within_sd

#### standardize with within SD ####
eb_climate_weekly_n <-
  eb_climate_weekly %>%
  melt(id.vars = id_vars, measure.vars = in_vars) %>%
  merge(within_sd, by = "variable") %>%
  .[,
    `:=`(value = (value - mean(value, na.rm = TRUE)) / within_sd),
    by = .(variable)
  ] %>%
  dplyr::select(-within_sd) %>%
  dcast(year + unit + week ~ ..., value.var = "value") %>%
  merge(eb_climate_weekly[, ..all_other_vars, with = FALSE], by = id_vars)

#### subset data to region-years with more than one Eurobarometer ####
eb_climate_outcome <- copy(eb_climate_weekly_n)
eb_climate_outcome <-
  eb_climate_outcome[
    !is.na(cc_eu_concern_w) & !is.na(emm_climate_reg_4w),
    n_eb_nutsyear := .N,
    by = .(nutsyear)
  ] %>%
  .[n_eb_nutsyear > 1,]

#### make centered matrices for mediation models ####
center_vars <- c(
  "cc_eu_concern_w",
  grep("tmean", names(eb_climate_outcome), value = TRUE),
  grep("emm.*_4w$", names(eb_climate_outcome), value = TRUE)
  )

eb_climate_outcome_c <- 
  feols(
    f_feols(yvar = "norm_dist", xvars = center_vars, fevars = m_fv),
    eb_climate_outcome, demeaned = TRUE) %>% 
  .$X_demeaned %>% 
  as_tibble()

#with google
center_vars2 <- c(
  center_vars,
  grep("google.*_4w$", names(eb_climate_weekly), value = TRUE)
)
eb_climate_outcome_g <- 
  feols(
    f_feols(yvar = "norm_dist", xvars = center_vars2, fevars = m_fv),
    eb_climate_outcome, demeaned = TRUE) %>% 
  .$X_demeaned %>% 
  as_tibble()

#### centered with spatio-temporal polynomials
center_vars <- c(
  "cc_eu_concern_w",
  grep("tmean", names(eb_climate_outcome), value = TRUE),
  grep("emm.*_4w$", names(eb_climate_outcome), value = TRUE),
  grep("_year_factor", names(eb_climate_outcome), value = TRUE)
  ) %>% 
  .[!grepl("2014", .)]

eb_climate_outcome_c_sp <- 
  feols(
    f_feols(yvar = "norm_dist", xvars = center_vars, fevars = c("unit", "year", "season", "event")),
    eb_climate_outcome, demeaned = TRUE) %>% 
  .$X_demeaned %>% 
  as_tibble()

#with google
center_vars2 <- c(
  center_vars,
  grep("google.*_4w$", names(eb_climate_weekly), value = TRUE)
)
eb_climate_outcome_g_sp <- 
  feols(
    f_feols(yvar = "norm_dist", xvars = center_vars2, fevars = c("unit", "year", "season", "event")),
    eb_climate_outcome, demeaned = TRUE) %>% 
  .$X_demeaned %>% 
  as_tibble()

#### save ####
save(
  eb_climate_weekly_n, # normalized panel
  eb_climate_outcome, # subset for outcome equations
  eb_climate_outcome_c, eb_climate_outcome_g, # centered data for use with lm()
  eb_climate_outcome_c_sp, eb_climate_outcome_g_sp, # centered data with polys for use with lm()
  file = "Data/eb_climate_reg.RData"
  )
